Aggregation of chemotactic organisms in a differential flow 
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We study the effect of advection on the aggregation and pattern formation in chemotactic systems 
described by Keller- Segel type models. The evolution of small perturbations is studied analytically 
in the linear regime complemented by numerical simulations. We show that a uniform differential 
flow can significantly alter the spatial structure and dynamics of the chemotactic system. The flow 
leads to the formation of anisotropic aggregates that move following the direction of the flow, even 
when the chemotactic organisms are not directly advected by the flow. Sufficiently strong advection 
can stop the aggregation and coarsening process that is then restricted to the direction perpendicular 
to the flow. 
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I. INTRODUCTION 

Directed motion of microorganisms and cells in re- 
sponse to chemical signals - chemotaxis - plays an im- 
portant role in a wide range of biological processes in- 
cluding migration of white blood cells, cancer invasion 
[ijl, embryonic development or in locating nutrients by 
bacteria, algae etc [H, [sl . In many cases the chemotactic 
cells not only detect, but also produce chemical signals 
that may attract other members of the population. This 
type of communication based on chemoattractant odors 
or pheromones can control group behavior, aggregation, 
swarming and collective decisions (quorum sensing) in 
bacterial colonies [4], slime mold [5| or insect popula- 
tions. Often the medium into which the chemical signal 
is released is not stationary but is a moving fluid (e.g. 
air or water) while the chemotactic cells or organisms 
are not transported by the flow as their motility is re- 
stricted to crawling on a solid surface. For example, mi- 
croorganisms may attach to surfaces developing biofilms 
0] found in natural environments or bioreactors and on 
a wide variety of surfaces, including living tissues, pip- 
ings, and industrial or medical devices. The interface 
between a surface and an aqueous medium such as water 
or blood provides an ideal environment for the develop- 
ment of microorganisms. The growth and structure of 
biofilm communities is a complex process regulated by 
the properties of the cell surface, diverse characteristics 
of the medium, type of substratum, and hydrodynamics 
of the aqueous medium. The influence of hydrodynamics 
on biofilm structures has been studied recently 0, 0, 
and was shown that the current velocity affects the struc- 
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ture and dynamics of natural biofilms resulting in differ- 
ent colony shapes Q. In Ref. [10] is described the use 
of "slow" laminar flows (from 100 /im s~^) to pattern 
cell culture substrate in capillary systems. Another in- 
teresting property which may influence the structure of 
biofilms is the cell-to-cell signaling or quorum sensing. 
For example, in Ref. |ll|] the importance of intercellular 
molecule signaling on biofilm differentiation is studied. 

The response of attached cells to a shear flow and the 
effects of cell-to-cell signaling on the aggregation have 
been also studied for a particular slime mold, the Dic- 
tyostelium discoideum. Using a laminar flow Decave and 
co-workers have established the critical shear stress for 
D. discoideum cells on glass [l2| and studied the mech- 
anisms responsible for the induced enhanced motility 
[l^. On the oder hand the aggregation of D. discoideum 
by means of secreting cyclic adenosine monophosphate 
(cAMP) has been modeled using stochastic and discrete 
approaches, the continuum descriptions of cell aggrega- 
tion have been mostly employed and later derived from 
mechanistic/microscopic descriptions [l^. The mathe- 
matical properties of these equations is relevant for a 
broad range models that have been developed to under- 
stand the aggregation process in a variety of organisms, 
pigmentation patterning, neural crest migration, inflam- 
matory response, tumor growth, etc. 

In this work we will study the simultaneous effect of 
differential advection and cell-to cell signaling on the ag- 
gregation and pattern formation of chemotactic biological 
populations using a model of partial differential equations 
to describe the evolution of the cell density and the chem- 
ical signal concentration. The resulting system is similar 
to the non-linear chemical reactions studied by Rovinsky 
and Menzinger involving activator and inhibitor kinetics 
where a differential flow can induce a pattern forming 
instability 
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II. MODEL 

A well known classical continuum model of chemotaxis 
at the population level is the Keller-Segel (KS) model [H, 
that describes the evolution of the density of chemotac- 
tic cells, i^(x, t), and the chemoattractant concentration, 
f(x, t), at point x and time t. When the chemical field 
V is advected by a uniform fiow V and the density field 
evolves on a fixed substrate we have 

dtu = \/' {DuVu - x{u)Vv) , (1) 
dtv = DyV'^v fu- sv - V - Vv. (2) 

where Du and Dy are constant diffusivities. The 
chemoattractant is assumed to be produced proportion- 
ally to the local cell density (with a constant of propor- 
tionality /) while it is degraded with a frequency s. Al- 
though in Ref. [I^ it was shown that a shear fiow in- 
creases the cell motility in D. discoideum, since we will 
consider slower fiow velocities than in Ref. [l^ (which 
were of the order of the detachment velocity), in Eq. [T] 
we assume that representative values of Du are closer to 
those measured in absence of fiow such as in Ref. [17|]. 
Assuming no-fiux or periodic boundary conditions the 
total mass of the biological component is conserved and 
can be characterized by the average density. In the orig- 
inal KS model [Sj, J^] V = and the chemotactic fiux is 
proportional to the particle density, i.e. x(^) = Xo^- Ex- 
tensions of this model with more general forms for x{^) 
have also been studied such as the chemotaxis model with 
prevention of overcrowding introduced in Ref. [l^ where 
X(^) Xou{l - u/umax) with Umax the maximum al- 
lowed cell density. An important feature of the KS model 
(observed in biological systems such as slime mold pop- 
ulations) is that it demonstrates an aggregation instabil- 
ity when the total mass of cells is larger then a certain 
threshold. Properties of the solutions of the KS system 
and its variants have been studied extensively (for recent 
reviews see Refs. 013, El). Interesting analogies be- 
tween KS type chemotaxis models and nonlinear mean 
field Fokker-Planck equations and generalized thermody- 
namics have been pointed out in [221 ]. 

In order to simplify the analysis of ([T]) and ([2j) we 
introduce non-dimensional variables by rescaling x' 
{s / DyY^'^x^ t' ^ 5t, u' ^ Uq^u^ and v' ^ s{fuo)~^v^ 
with uo the initial mean cell density, resulting the follow- 
ing system: 

dtu = \/ ■ {DVu - x(^)Vv) , (3) 
dtV = V'^v -^u-v-W-Vv, (4) 

where we have defined D = Du/Dy^ xi'^Y = 
f /{sDy)x{u), V = {sDy)~^^'^\ and omitted primes. In 
order to estimate the typical spatial and temporal scales 
in this problem we can use parameter values given in Ref. 
[Tst for the chemoattractant cAMP: Dy 10"^ cm^ s"^ 
and s ~ 1 s. Thus the typical units for the rescaled 
length, time and velocity are of the order of 0.1 /im, 1 s 
and 10 /im s~^ respectively. 



It is important to note that for D. discoideum the crit- 
ical shear stress for detachment on glass is of the order of 
(7i/2 ^ 2.6 Pa [12]. For low Reynolds number when the 
inertial effects can be neglected the wall shear stress on 
the adhering cells is proportional to the uniform velocity 
following (J = 6r]V/d^ where rj is the dynamic viscosity 
of the fiuid and d is the distance between the top of the 
chamber and the substrate. Using d = 0.25 mm as in 
Ref. [13 and water at room temperature {r] = 10~^ Pa 
s) we obtain an estimate for the detachment velocity: 
Vi/2 ^ 10^ Thus, as we will see below, the fiow 

velocities considered here are smaller than the velocity 
needed to detach the cells from the substrate and we can 
assume that the organisms are not advected by the fiow, 
although they are still able to move by crawling on the 
solid surface as represented by the chemotactic and dif- 
fusive terms. 



III. LINEAR ANALYSIS 

In order to gain insight into the system we consider 
the stability of the spatially uniform solution. Assuming 
uniform initial conditions for u and v with v{t = 0) = vq 
all spatial derivatives vanish in Eqs. (|3]) and (|4j) and we 
have the solutions u{t) = 1 and v{t) = 1 -\- {vq — l)e~*. 
Thus, independently of the initial conditions the con- 
centration tends to 1 for large times. To investigate 
pattern forming instability in this system we consider 
the evolution of spatially non-uniform periodic pertur- 
bations added to the uniform steady state of the form 
ii(x, t) = 1 + 'uexp [iq • X + cc;(q)t] and '^(x, t) = 1 + 
exp [iq • X + C(;(q)t], and study whether the perturba- 
tion is amplified or damped out in the course of time. 
Here q is the wave vector of the perturbation, Cc;(q) is 
the corresponding dispersion relation, and u and v are 
the amplitudes of the perturbation at t = 0. Substi- 
tuting these expressions into (|3]) and (|4]), and neglecting 
quadratic terms in we obtain a linear system of 

equations where non-trivial solutions only exist if the de- 
terminant of the coefiicient matrix is equal to zero. In 
contrast to previous chemotactic models where advection 
is not considered (see for example [23| for a linear stabil- 
ity study of an inertial model generalizing the KS model), 
the resulting quadratic equation for the dispersion rela- 
tion has complex coefficients and reads 

cj^ + uj{a + ib) + (c + id) = 0. (5) 

The coefficients a, 6, c, and d are functions of parameters 
and wave- vector components yielding 

a = (L> + l)q^ + 1 (6) 
6 = V • q (7) 

c = Dq^^[D-x{l)]q' (8) 
d = Dq^\ ' q. (9) 
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The real and imaginary parts corresponding to the two 
complex solutions are (see for example page 95 of 

Re(c.±) = -l± ^{[(a^ -h'- 4c)^ + (2ah - Adf]'^^ 

-^a'-h'-AcY^^ (10) 

T / ±\ ^1 sgn(2a6-4(i) ^ i2 a \2 
Im(cj±) = --± ^ ^{[(a^ -b^- 4c)2 

+ {2ab - 4d)2] V' - + 6^ + 4cy/\ (11) 

The real part of co gives the growth or decay rate of the 
perturbation amplitude. In particular, the mode corre- 
sponding to the maximum of Re (a;), which we denote 
by q^, determines the characteristic wavelength of the 
pattern in the linear regime, while the imaginary part 
describes its propagation in space. The velocity at 
which the instability travels across the substrate, corre- 
sponding to the phase velocity of the mode q^, satisfies 
the relationship [25] 



. q^ = -Im [u;{q[ 



(12) 



The negative branch of the dispersion relation is un- 
conditionally stable, i.e. Re(cc;~) < 0, but the positive 
branch may produce non-trivial dynamics and pattern 
formation for a certain range of wavenumbers. Some in- 
sight into the behavior of the system can be obtained by 
investigating the limiting cases of small and large wave- 
lengths. For large values of q we can expand Eq. ([TQ|) 
to obtain Re(cc;+) = —Dq^^ therefore the amplitude of 
perturbations decays exponentially for these modes. In 
the case of large wavelength perturbations the expansion 
to the lowest orders in q yields 



Re(a;+) = [x(l) - D] 
- X{l)q' {[I - D + x{l)] 



(V-qf} + 0(q6). 

(13) 



Thus, when x(l) — ^ > the positive branch has a band 
of unstable modes for small wave vectors. In the case of 
the standard KS model, x(l) = Xo, the above condition 
is equivalent to the aggregation threshold: D/xo < 1 
0, In Eq. ^ we observe that the advection ve- 

locity appears at the order q^. This means that, when 
the above condition for instability is satisfied, there is 
always a band of unstable modes with long wavelengths 
around <7 = 0, but the range of unstable modes and their 
growth rate decreases with the advection velocity. This 
stabilizing effect of the fiow is in contrast with the be- 
havior of reaction-diffusion systems studied by Rovinsky 
and Menzinger [l5|, where the differential fiow induces 
an instability at a finite wavelength. For the imaginary 
part of cj, we have the following expansion of Eq. ([TT]) for 
small wave vectors: 



lm{u;+) = -x{l)q^ (V • q) + O (g^ 



(14) 



chemotaxis induces a phase velocity, V^, that is inversely 
proportional to the square of the wavelength and, conse- 
quently, it is small relative to the advection velocity. 



IV. PATTERN FORMATION IN ONE 
DIMENSION 

Figure [Tfa) shows the effect of the advection velocity 
on the real part of the dispersion relation given by Eq. 
(p!Q|) for a ID system. When V increases the wavenum- 
ber of the dominant mode corresponding to the maxi- 
mum of Re(cc;+), q\ decreases [Fig. [Hb)]. In fact, from 
the expansion of Eq. (fTQ|) for small q given by Eq. (p!3|) 
we have q^ = [u/ (2/C)]"^^^ with u = x(l) — D and 
JC = x(l) [(1 - ^) + X(l) + V^] . Thus for large V, the 
dominant wavelength in the linear regime, = 27r/q\ is 
proportional to the advection velocity. The phase veloc- 
ity of the spatial pattern [Eq. ([12])] is also shown in Fig. 
[TJb) as a function of the advection velocity for the exact 
dispersion relation [Eq. (pT]) ]. We notice that the phase 
velocity has a maximum for a certain value of V. Thus, 
surprisingly, when the advection velocity is increased be- 
yond this value, the phase velocity of the pattern de- 
creases as V is increased. This non-monotonic depen- 
dence can also be shown using the expansion of lm{u;~^) 
for small wavenumber modes [Eq. ([T4j) ] from where a 
compact analytical expression for the phase velocity is 
obtained [see Fig. [Hb)] 



v' = x{i)v{qy 



[X(l) -D]V 



2[{l-D)+x{l) + V^Y 



(15) 



Therefore, although the particle density is not directly 
advected by the fiow, using Eq. (fT2]) we see that the 



Although this expression is only valid for small values of 
q\ it is qualitatively similar to the exact solution and 
already shows that increases linearly for small values 
of V and when th e advection vel ocity exceeds a certain 
threshold, Vt = ^ {1 — D) -\- x{^)j larger values of V in- 
duces slower pattern movement. 

It is interesting to discuss how the linear wavelength 
and the phase velocity are modified when D is decreased, 
since in typical experiments Du/Dy <C 1. For strong ad- 
vection we can use the expansion of Eq. (p!Q|) for small 
q to obtain the minimal value of the linear wavelength. 
Since in the aggregation regime, < D < x{^)^ the lin- 
ear wavelength is an increasing function of for D ^ 
we obtain A^^^ = + x(l) + from where we 

see that the linear wavelength becomes larger when the 
chemosensitivity is increased. On the other hand, Vt is a 
decreasing function of D. Thus, in the limit D ^ ^ the 
maximum phase ve locity is reached when the advection 
is = ^JY^^^xSX)• Therefore, as we could expect, the 
efficiency of the particles following the chemical field de- 
pends on the chemosensitivity. For larger values of x(l) 
the cells can more accurately follow the chemical field up 
to larger values of the advection velocity. 

For the numerical simulations we use the chemotactic 
response function: x('^) = Xo'^(l — u/umax)i that avoids 
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FIG. 1: (a) Real part of the dispersion relation Re(a;^) given 
by Eq. (|1Q|) as a function of q for a ID system with x{^) — 
Xou{l-u/umax), for D = 1, Xo = 2.5, Umax = 4 and different 
values of the advection velocity, V. (b) (black lines) and 
(red Unes) as a function of the advection velocity. The 
solid lines represent the exact solution obtained from ([T0|) and 
(|lip . The dashed lines represent the approximate solution 
for small wavelengths. The velocity of the pattern and the 
dominant mode measured from the numerical simulations are 
represented by triangles and circles respectively. 



the singularities associated with other models. Interest- 
ingly, we have also found that in the case of the standard 
linear response, x('^) = Xo'^, advection can suppress the 
singularity and produce qualitatively similar behavior to 
the previous function when V is sufficiently large. We 
consider periodic boundary conditions with initial fields 
in the uniform steady state with a small amplitude ran- 
dom perturbation. The parameters in Eqs. (|3]) and (|4|) 
were set to D = 1, xo = 2.5, and Umax = 4 for the one- 
and two-dimensional simulations. 

Representative one-dimensional numerical simulations 
of the particle density profiles u{x^ t) are shown in Fig. [2] 
for two different values of the advection velocity. In the 
case oi V = 1 [Fig. Efa)] we observe a slow coarsening 
process which tends to reduce the number of structures as 
observed in the system without advection However, 
for V = 5 [Fig. EJb)] after a short transient a periodic 
pattern with constant wavelength develops. As shown in 
Fig. [21 in the presence of advection the x —x symme- 
try of the profiles is broken [see also Fig. [31(a)] and the 
pattern propagates in the positive direction (i.e. in the 
direction of the advection velocity). In the absence of 
coarsening (i.e. for values of F > 2) the pattern moves 
uniformly with a well defined velocity as shown in Fig. 
[Sfb). The velocity of the pattern was measured for differ- 
ent values of V and is plotted in Fig.lljb). The measured 
velocity agrees well with the analytical results obtained 
from the linear stability analysis. Independently of the 
initial chemical concentration, after a transient time the 
behavior of the chemoattractant v is very similar to the 
density profiles, as shown in Fig. [31(a) where the pro- 
files of u and v are plotted together for different values 
of V. For large advection velocities the particles can 
not "follow" the chemical gradients, the two profiles are 
more different and the aggregates are more spread out. 
Figure [S^b) shows the temporal evolution of the wave- 



length X{t) (measured as two times the average distance 
between consecutive minima and maxima) for different 
values of V. Aggregation starts earlier for larger values 
of V and for strong enough advection values (larger than 
Vt) coarsening is interrupted and the wavelength of the 
pattern becomes constant with a final value which is pro- 
portional to V as predicted by the linear analysis [see Fig. 
mb)]. Thus the nonlinear effects (such as coarsening) are 
avoided for large values of V and the system remains in 
the linear regime. Similar results are obtained for the 
standard KS model, x('^) = Xo'^? with the difference that 
for slow advection the numerical simulations do not reach 
a stationary state, indicating an aggregation singularity 
with unlimited growth of the density in some points. 





FIG. 2: (a) Particle density profiles at equally spaced times 
(profiles are offset vertically) between t = (bottom) and 
t = 3500 (top) for: (a) V = 1 and (b) F = 5. 





FIG. 3: (a) Density, (solid line) and signal concentration, 
V, (dashed line) profiles at t = 3500 for V — 0,1,2,5, and 
10, from bottom to top respectively (profiles are offset ver- 
tically), (b) Temporal evolution of the lateral pattern wave- 
length, A(t), for different values of V. Error bars represent 
the standard deviation calculated from 150 different random 
initial conditions. 



V. PATTERN FORMATION IN TWO 
DIMENSIONS 

The effect of the advection term in the real part of the 
dispersion relation given by Eq. ([10]) for a 2D system is 
shown in Fig. [H Without loss of generality we assume 
that the advection is along the x-axis. Thus, we can write 
V = (F, 0) and we plot Re(cc;+) as a function of qx and 
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qy for different values of V . When ^ = the system 
is isotropic and ah the wave vectors within a distance 
from the origin maximize Re(c<;+) [Fig. [41(a)]. When V is 
nonzero, the maximum is oriented along the ^-axis but 
the absolute value is not affected by the value of V [Fig. 
[11(b)]. The effect of the advection on the wavelength and 
orientation of the pattern can be analyzed in the small q 
limit using Eq. ([T3|) . It can be shown that the dominant 
mode is always oriented along the ^-axis and its value is 
ql = [u/ where /C, = x(l) [(1 - D) + x(l)] and 

TJ is the same as defined for the ID case (see Appendix). 
Thus, in contrast to the one-dimensional system, the ob- 
served dominant wavelength of the pattern in the linear 
regime is independent of the advection velocity V ^ al- 
though the features of the pattern are altered by the 
symmetry breaking induced by the advection. For the 
phase velocity, using the expansion of Im(a;+) for small 
q and Eq. ^ we have • q^/q^ = x(l)g^ (V • q^) = 
because V = (y,0) and = (0,g^). Therefore, the 
dominant mode in the pattern is not moving since it is 
oriented in the y axis that is perpendicular to the advec- 
tion velocity. However, for modes in the x-axis, qx, the 
phase velocity in this direction is = x(l)^^^ which is 
proportional to the advection velocity and increases with 
the wavenumber. 




FIG. 4: Real part of the dispersion relation Re(a;^) given by 
Eq. (|lQp as a function of Qx and Qy for a 2D system with (a) 
y = and (b) V = 2. 



Two-dimensional numerical simulations were per- 
formed using a cubic interpolation semi-lagrangian 
scheme [26^ for the advection. For the spatial discretiza- 
tion of the particle density we used the method proposed 
recently by Grima and Newman [23], that allows for an 
accurate analysis of the evolution of the system without 
the dissipative effects of other schemes. In the absence 
of advection a slow coarsening process occurs in which, 
after long times, the organisms accumulate into a single 
aggregate [l9| . When advection takes place the resulting 
pattern is not isotropic and the dominant wave vector is 
oriented along the y axis. The temporal evolution of u 
for different values of the advection velocity is shown in 
Fig. [5l As in the ID case, the evolution of the chemical 
field is very similar to the density field (see supplemen- 



tary movies (Hi)- As predicted above, the pattern does 
not propagate in the y direction and the coarsening pro- 
cess along the y axis is not affected by the advection 
(see supplementary movies [HI). Thus, in this direction 
the pattern features remain the same when the advec- 
tion velocity is increased. A different scenario is found in 
the X direction where the pattern moves with a velocity 
which increases with V. The simulations also show that 
smaller aggregates move faster than larger ones, consis- 
tently with the dependence of the phase velocity on the 
wave-number in the linear analysis. As in the ID case, 
the characteristic wavelength of the pattern in the x di- 
rection increases with V. For slow advection there is a 
clear coarsening which leads to smaller propagation ve- 
locity of the aggregate as predicted above. For finite 
system sizes and large values of the advection velocity, 
the wavelength of the dominant linear instability in the 
X direction becomes larger than the system size and the 
resulting pattern is homogeneous in the x direction, ob- 
serving a static pattern in this direction (see Fig. [5l for 
V = 10). 




FIG. 5: Temporal evolution of u for different values of the 
advection velocity: V = 1 (first row), V = 2 (second row), 
y = 10 (third row); and different times: t = 35 (first column), 
t = 70 (second column), t = 200 (third column), t = 870 
(fourth column). The x axis is horizontally oriented and the 
system size is 128 x 128. See also supplementary movies 



VI. CONCLUSIONS 

In this work we have studied theoretically how an ad- 
vected uniform flow influences the aggregation dynam- 
ics in Keller-Segel type models. We found that, in the 
presence of a differential flow, an advective instability 
produces a pattern moving in the direction of the flow. 
Interestingly, although the organisms are not directly ad- 
vected, the advection of the chemotactic signal induces a 
movement of the particle density as the organisms try to 
follow regions of high chemical concentrations. When the 
organisms mobility due to chemotaxis is weak in compari- 
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son to the advective flow (which is typicahy much faster), 
the balance between the chemotactic flux of the organ- 
ism density and the advective transport of the chemical 
fleld breaks down and it is then restored by a change in 
the characteristic of the spatial patterns, that become 
strongly anisotropic with elongated stripe-like structures 
aligned to the direction of the flow. Furthermore, as 
shown above, for flows larger than a certain threshold 
the organism can not follow the chemotaxis signal reduc- 
ing their velocity and eventually preventing the forma- 
tion of large aggregates. Thus, the presence of a linear 
advection term inhibits the formation of large gradients 
and diminish the nonlinear effects. This stabilizing effect 
may completely stop the chemotactic aggregation and 
coarsening process in the direction of the flow. Although 
a simple uni-directional flow can not suppress the coars- 
ening in the perpendicular direction, preliminary results 
with more general non-uniform time-dependent velocity 
flelds in 2D show that aggregation may be halted pre- 
venting the appearance of singularities associated to the 
KS models. This is similar to the arrested coarsening 
process observed in binary mixtures [1^, [sO] • 

The theoretical results on the distribution of chemotac- 
tic cell populations under the influence of a differential 
flow are relevant for various natural and artiflcial systems 
including biofllms and could also be studied experimen- 
tally in the context of Dictyostelium aggregation. This 
work may also provide a starting point for the study of 
more general biological pattern formation phenomena in 
advected environments. 



value of the real part of the dispersion relation and its 
wavelength is associated to the wave vector = (ql^Qy). 
This vector verifles 



(Al) 







■5Re(cc;+)' 


dqx 







which have the following real independent solutions 



qo = (0, 0) , qi = (y^, 0) , q2 = [q, ^) ,(A2) 



where z^, /C, and Ky are deflned as in the main text and 
assumed positive. In order to decide which of the remain- 
ing solutions provide the absolute maximum of Re(cc;+), 
we finally substitute the wave vectors given by (|A2p into 
Eq. ([T3|) : we obtain simply 



[Re(a.+)]^^ = 0, = ^, [Re(^+)] 



'i^ 4/C,' 



(A3) 

from where, since /C > /C^y ^ 7^ 0^ we conclude that 



q2. 
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